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Neutrinos in core collapse supernovae are likely trapped by neutrino-nucleus elastic scattering. 
Using molecular dynamics simulations, we calculate neutrino mean free paths and ion-ion correla- 
tion functions for heterogeneous plasmas. Mean free paths are systematically shorter in plasmas 
containing a mixture of ions compared to a plasma composed of a single ion species. This is be- 
cause neutrinos can scatter from concentration fluctuations. The dynamical response function of 
a heterogeneous plasma is found to have an extra peak at low energies describing the diffusion of 
concentration fluctuations. Our exact molecular dynamics results for the static structure factor 
reduce to the Debye Huckel approximation, but only in the limit of very low momentum transfers. 
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I. INTRODUCTION 

Core collapse supernovae are extraordinarily energetic 
explosions where 99% of the energy released is radiated 
in neutrinos P, This is because weakly interacting 
neutrinos are the only known particles that can diffuse 
quickly out of the dense stellar core. When the density of 
the core reaches w 10 12 g/cm 3 , neutrino nucleus elastic 
scattering is thought to temporarily trap the neutrinos 
3]. The neutrinos then provide a degeneracy pressure 
that helps support the star. If the neutrinos were not 
trapped, the star could collapse to a black hole without 
a large supernova explosion. Thus, the dynamics of a 
supernova is sensitive to how neutrinos interact in dense 
matter. For example, changes in the electron capture 
rate on nuclei could change the sensitive balance between 
electron Fermi pressure and gravity 0, Q . However elec- 
tron capture can only proceed if the produced neutrinos 
are not Pauli blocked. 

The neutrino wave length is comparable to the spacing 
between ions. Therefore, ion-ion correlations can signifi- 
cantly modify neutrino-nucleus scattering cross sections. 
An ion in a plasma will be surrounded by a screening 
cloud of other ions and perhaps electrons. Neutrino inter- 
actions with this cloud will, in general, reduce neutrino- 
nucleus cross sections and increase the neutrino mean 
free path, see for example . Most calculations of this 
screening assume a one component plasma composed of a 
single species of ion, see for example f7| . A first attempt 
was made in ref. 8] to consider scattering from a plasma 
including both alpha particles and heavy nuclei. How- 
ever, this was based on a simple prescription to describe 
the mixture. 
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If the plasma has a mixture of species with different ra- 
tios of weak to electromagnetic charges then the screen- 
ing of neutrino reactions can be very different. Neutri- 
nos scatter from fluctuations in the weak charge den- 
sity. Fluctuations in composition, that change the weak 
charge density but do not change the (electric) charge 
density, will feel no electric restoring force. Therefore, 
these fluctuations may not be screened and can lead to a 
large change in the neutrino cross section. This was first 
discussed in connection with ion and electron screening 
in refs. Recently Sawyer argued that com- 

position fluctuations for a mixture of ions with different 
ratios of neutron number N to charge Z can increase the 
cross section. However Sawyer used the Debye Huckel ap- 
proximation which is only valid at very low momentum 
transfers. 

In this paper we use molecular dynamics simulations 
to accurately calculate neutrino scattering cross sections 
and mean free paths from a dense plasma composed of 
mixtures of different ions. Section II discusses the mix- 
ture of ions expected in supernova plasmas and presents 
our neutrino scattering formalism, section III describes 
our simulations, section IV our results and we conclude 
in section V. 



II. FORMALISM 

We are interested in neutrino scattering from a dense 
plasma during the infall phase of a supernova. As the 
density and temperature rise the rate of nuclear reac- 
tions increases dramatically. This should bring the com- 
position of the plasma into nuclear statistical equilibrium 
(NSE). In NSE all nuclei are in chemical equilibrium 
so the abundance of a given state is determined by its 
binding energy and entropy. We first discuss this NSE 
composition and then describe how to calculate neutrino 
scattering from this heterogeneous mixture. 
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A. Composition 

The composition of a plasma in nuclear statistical equi- 
librium depends on the density, temperature and proton 
fraction. For given conditions, we expect the most abun- 
dant heavy isotope to have a binding energy and entropy 
that leads to a large free energy. However nuclear reac- 
tions can quickly add or subtract nucleons to this most 
abundant species. Therefore, we expect a distribution of 
isotopes near the most abundant one. In addition, en- 
tropy can favor some light isotopes such as alpha parti- 
cles and free neutrons. There have been many statistical 
models that calculate the NSE composition, see for ex- 
ample 12]. Typically, these statistical approaches incor- 
porate detailed binding energy models but only include 
the interactions between nuclei approximately. 

Alternativel y, a semi-classical microscopic model de- 
scribed in ref. [l3Lll4| includes the strong interactions be- 
tween nucleons in different nuclei in an identical fashion 
to the interactions between nucleons in a given nucleus. 
This allows the model to describe not only a plasma of 
isolated ions but higher density matter where the nuclei 
merge together to form complex pasta shapes and then 
uniform nuclear matter. 

However, this simple model may not reproduce the 
binding energy of individual nuclei as well as statisti- 
cal models. Although the force has been fit to reproduce 
the binding energy of nuclear matter, the model does not 
include pairing and shell structure. Comparing the mi- 
croscopic and statistical models provides some estimate 
of the possible range in compositions that could be ex- 
pected. 

In the microscopic model, the location of each nucleon 
is followed in a molecular dynamics simulation. For ex- 
ample, in ref. Ti} trajectories for 40,000 nucleons have 
been calculated by integrating Newtons laws. At any in- 
stant in time, the configuration of nucleons is divided into 
nuclei using the following simple algorithm. A nucleon is 
said to belong to a given nucleus if it is within a cutoff 
radius R cu t = 3 fm of at least one other nucleon in the 
nucleus. This algorithm uniquely divides a given config- 
uration of many nucleons into a collection of nuclei. The 
resulting distribution of nuclei will be presented in sec- 
tion llVI However, first we present our neutrino scattering 
formalism. 



B. Neutrino Scattering 

In this section we describe how neutrino scattering 
is modified by ion-ion correlations. The free neutrino- 
nucleus elastic scattering cross section is 

d g /m= G2c2E ^ 2 +cose \ (1) 

Here G is the Fermi constant, E v is the neutrino energy, 
the scattering angle and C is the total weak charge of 



a nucleus with charge Z and neutron number N, 

C = -2Z sin 2 O w + (Z - N)/2, (2) 

with a Weinberg angle of sin 2 <dw = 0.223. Note that 
sin 2 <dw ~ 0.25. In the following we approximate C « 
—N/2. For a mixture of ions we will use, 

C = -\{N), (3) 
where the average neutron number (N) is, 

( N ) = N-H ( 4 ) 

for a system of Ni on total ions where the ith ion has 
neutron number Ni. 

Ion correlations can be taken into account by multi- 
plying dao/dQ by the static structure factor S(q) [l5| . 

da/dn = da /dQS{q). (5) 

Here q is the momentum transfer and <jq is the free cross 
section. The static structure factor adds coherently the 
contributions for neutrino scattering from different nu- 
clei, including the relative phases, and can be calculated 
from, 

S(q) = -^((pH<im))-\(p(<l))\ 2 ), (6) 
with /5(q) the density operator is given by, 

/5(q) = ex p(*q ■ r *)- ( 7 ) 

Note the choice of the normalization (N) in Eqs. Ij3l7|) is 
a somewhat arbitrary convention. We make this choice 
because one often approximates a mixed system with a 
single species where all the Ni = (N). Alternatively one 
could replace (N) by (N 2 ) 1/2 in both Eqs. © and 
with no change in the cross section in Eq. J3J). With this 
second choice S(q) — » 1 at high q. 

The transport mean-free path A is inversely propor- 
tional to the transport cross section <r*, A = 1/piCr* with 
Pi the ion density. In a medium, the transport cross sec- 
tion cr* can be obtained by multiplying the free transport 
cross section 

(Jq = J dQ(l - cos9)da Q /dn = ^G 2 C 2 El/ir (8) 
by (S), 

a* = *t(S), (9) 
where (5) is the angle average of S(q) 0, 



(S{E v ,p,T)) = 




dcos (9(1 + cos 0) (1 - cos 8)S{q{9)). 



(10) 

Here q{9) 2 = 2E 2 V {\ - cos0) and the factor of (1 + cos0) 
comes from angular dependence of the free cross section, 
Eq. (JTJ. 
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C. Debye Huckel Approximation 



III. SIMULATIONS 



At very low momentum transfers q the simple Debye 
Huckel approximation is valid [Ti| . This provides insight 
into our full molecular dynamics simulation results. The 
static structure factor in the Debye Huckel approxima- 
tion Sg H is given by 



nDH _ ±_ 

q Mm* 



i,j=l 

Here T is the temperature and 

13 - ij q 2 + q 2 + e 2 £ m „ Z m Z n Tl mn ' 



(11) 



(12) 



where electron screening is describe by q 2 = (4a/ir)kp. 
In the Debye Huckel approximation, 



tt _ 3jj n i _ 

LHj - T - y T , 



(13) 



where we assume the density of the ith ion is n$ = l/V 
with V the system volume. The total ion density is pi — 
Ni on /V . Using this expression for reduces Kij to, 



y VT \ y yT q2 



K 2 J 



with K = K~ on + q e and 



The final result is 



VT 



q (N) 2 



zZ z * 



where, 



(A fc ) 



Z.z? 



for k — 1,2 and, 



1 » i nn. 



(14) 



(15) 



(16) 



(17) 



(18) 



In the <7 = limit, and neglecting electron screening q 2 — 
0, we have 



(19) 



Thus the static structure factor depends on the disper- 
sion in the ratio of weak to electromagnetic charge, with 
(AA) 2 = (A 2 ) — (A) 2 . Alternatively for a single compo- 
nent plasma, we have 



S. 



dh _q 2 + ql 



q 2 + k 



(20) 



We compare our full simulation results to the Debye 
Huckel expression, Eq. (|16|l . in Section llVl 



In order to calculate the static struture factor, for ar- 
bitrary momentum transfer q, we perform molecular dy- 
namics simulations using the Verlet algorithm ^tJ- F° r 
the conditions we consider, the thermal deBrogile wave- 
length of the ions is much less than the inter-ion spac- 
ing so we assume that the ions behave classically. We 
are interested in momentum transfers much less than the 
electron Fermi momentum q << kp. Therefore, we can 
describe the interaction between ions with a Yukawa po- 
tential El: 



V(i,j) 



(21) 



where is the distance between a pair of ions. The 
electron screening length is A e = tt/ (ekp) where the elec- 
tron Fermi momentum is kF = (peS^ 2 ) 1 ^ 3 , the electron 
density is p e — (Z)pi, and e 2 = Aira. 

Periodic boundary conditions are used to minimize fi- 
nite size effects. The distance between ions is then 
given in terms of the coordinates x,y and z of the ith 
and jth particles in the form: 



Tij = y/lxi-xtf + lvi-ytf + lzi-Zj] 2 . (22) 

Then, the distance used for the periodic boundary con- 
dition is: 



min(|Z|,Z- \l\). 



(23) 



Here L is the side of the cubic box V containing the JVj on 
ions, L = V 1 / 3 . For ion density pi the box volume is 
V = Nion/pi. 

We used three different Fortran molecular dynamics 
programs to do our simulations. The programs are indi- 
cated by the author's initials. All the codes use the direct 
particle-particle method to calcluate interactions, where 
one directly sums the ^N(N — 1) interactions among N 
ions. Thus the amount of work to do a simulation in- 
creases as the square of the number of particles. 

The CH and LC codes are serial codes suitable for do- 
ing small problems. They were used for MD simulations 
of 1000 to 4000 ions. The larger 10000 and 40000 ion runs 
however were too much work for a serial code, so these 
were done with a parallel program, DB. The 40000 ion 
mixture run was particularly time-consuming, both be- 
cause of the 100 fold increase in work per time step over 
the 4000 ion run, and because it had to be run longer 
in order to capture the diffusion of the weak charge fluc- 
tuation across the larger simulation box. DB uses the 
MPI (Message Passing Interface) library to pass messages 
among processes. The iV ions are partitioned into p sets, 
where p is the number of MPI processes, and each pro- 
cess is tasked with calculating the forces on its set of N/p 
particles due to all N. This makes for ^N(N-l) jp inter- 
actions per process per time step. The main communica- 
tion needed is for all processes to share their coordinates 
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with each other, what is called an "allgather" . This can 
result in a high communication overhead, which is rel- 
atively less costly for large problem sizes. After swap- 
ping coordinates, the force calculation and integration of 
Newton's equations can proceed without communication. 
The force calculation in each MPI process is addition- 
ally parallelized with OpenMP, which is a set of compiler 
directives used primarily for parallelizing DO loops on 
multiprocessor shared memory machines. 

Run DB1 was done on a 36 node distributed memory 
parallel computer at Indiana University, the AVIDD-0 
machine, where each node board has two AMD Opteron 
processors (Advanced Micro Devices). Each node ran 
one MPI process consisting of two OpenMP threads. We 
started with the 40000 ions uniformly and randomly dis- 
tributed in a simulation box of edge length L = 822.8 fm, 
and with velocities distributed according to a Boltzmann 
distribution at temperature T = 1.00 MeV. We brought 
the system to equilibrium during the first 7000 warmup 
steps, with a time step At = 25 fm/c, and then ran a fur- 
ther 252500 measurement steps, with At = 50 fm/c. Ev- 
ery tenth configuration was written to disk. These were 
used for calculating S(q, w) and S(q). The velocity- Verlet 
algorithm was used to integrate Newton's equations. As 
this algorithm conserves total energy, we rescaled veloci- 
ties every 1250 time steps in order to maintain the tem- 
perature at T — l.OOMeV. We used 4 Opteron nodes (8 
processors) for the warmup steps, 8 nodes (16 processors) 
for the first 65000 measurement time steps, and 10 nodes 
(20 processors) for the remaining 187500. The program 
did not have perfect linear speedup with the number of 
nodes. A 40000 particle problem is actually too small to 
run efficiently on this number of processors, due to com- 
munication overhead. Speedup becomes more linear as 
the problem size is increased. The full DB1 run consisted 
of 259500 MD time steps, and took 793 hours spread out 
over 8 weeks. 

We did the 10000 pure ion run DB2 on one dual pro- 
cessor PowerPC970 node in Indiana University's 42 node 
IBM JS20 system. For this run we compiled the code 
to use only the OpenMP capability, as the 2 proces- 
sors share the board memory. We again started with 
the ions uniformly and randomly distributed in a simu- 
lation box of edge length L = 518.4 fm, and with veloci- 
ties distributed according to a Boltzmann distribution at 
T = 1.00 MeV. We did 8000 warmup steps at At = 25 
fm/c, and 24000 more at At = 50 fm/c. We then did 
100000 measurement steps at At — 50 fm/c, again writ- 
ing out every tenth configuration for analysis. Velocities 
were rescaled every 1250 time steps to maintain the tem- 
perature at T = 1.00 MeV. The full DB2 run consisted 
of 132000 MD time steps and took 346 hours spread out 
over 16 days. We did run DB2 almost simultaneously 
with DB1, thus running on two processors caused no ex- 
tra delay. 

Finally, we measure the static structure factor S(q) for 



valuse of q, 

{q x ,q v ,qz} = -j-(n x ,n y ,n z ) . (24) 

To minimize finite size effects we chose n x , n y , n z , as 
integers. We average over directions of q to improve the 
statistics. Statistical error bars for codes CH and DB are 
estimated from the dependence of S(q) on the direction 
of q. Statistical error bars for code LC are estimated by 
dividing the measurement time into a number of separate 
groups and looking at the variance of S(q) for different 
groups. 

IV. RESULTS 

In this section we present results for the static struc- 
ture factor from full molecular dynamics simulations 
and compare to the simple Debye Huckel approximation. 
However first we present results for the plasma composi- 
tion based either on a simple microscopic model or on a 
statistical model. 



A. Composition 

In this subsection we present results for the distribu- 
tion of charge Zj, and neutron number Ni to be used in 
our molecular dynamics simulations. We start with con- 
ditions considered in ref. ^4[: a density of 0.01 nucleons 
per fm 3, or 1.66 x 10 13 g/cm 3 , a temperature of T = 1 
MeV and a proton fraction Y p = 0.2. Reference 0] 
performed a molecular dynamics simulation in the nu- 
cleon coordinates with a simple short range nuclear plus 
Coulomb interaction. The nuclear interaction was fit to 
reproduce the saturation density and binding energy of 
nuclear matter. 

The final coordinates of a simulation with 40000 nucle- 
ons are divided into nuclei by assuming a nucleon belongs 
to a nucleus if it is within 3 fm of at least one other nu- 
cleon in the nucleus. This algorithm produces a collection 
of about 300 medium mass nuclei with (A) m 100 shown 
in Fig. ^ In addition there are a few light nuclei with 
mass A < 10 (not shown) and about 10000 free neutrons. 

We are interested in the effects of ion-ion correlations 
so we neglect the free neutrons. Neutrino interactions 
with these free neutrons are not expected to have large 
correlation effects. In addition, we also neglect the light 
nuclei. The light nuclei have small weak and electromag- 
netic charges, and are expected to play only a small role 
in neutrino scattering. Furthermore, the light nuclei have 
larger thermal velocities that require a smaller molecu- 
lar dynamics time step which slows down the molecular 
dynamics simulations. 

We compare these microscopic results to the statistical 
model of Botvina et al. 01 which we also show in Fig. 
^ For a slightly lower density 10 13 g/cm 3 and the same 
Y p and T, the statistical model yields significantly larger 
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TABLE I: Results for ion density pi, average mass number 
(A), average charge (Z), and mass fraction of free neutrons 
X n for microscopic [l4| and statistical fl^l models. 
Model p (gm/cm 3 ) Y p p r (fin" 3 ) (A) (Z) X n 

Microscopic 1.66 x 10 13 0.2 7.18 x 1(T 5 105.2 32.1 0.25 
Statistical 10 13 0.2 1.60 x 10" 5 328.6 75.7 0.33 

Statistical 10 12 0.5 1.06 x 10~ 5 56.8 28.3 « 



In a core collapse supernova the proton fraction starts 
out just below Y p — 0.5 and decreases with time due 
to electron capture. Therefore, in Figs. |3 El we show 
statistical model results for Y p = 1/2, T = 1 MeV and a 
density of 10 12 g/cm 3 . Note that the microscopic model 
has not been simulated for Y p = 1/2. The dispersion in 
N/Z is smaller for Y p = 1/2 than for Y p = 0.2. 



FIG. 1: Distribution of charge Z versus neutron number N 
for ions at a temperature of T = 1 MeV and proton fraction 
Yp — 0.2. The microscopic Pasta results are from ref. |l4| at 
a density of 1.66 x 10 13 g/cm 3 while statistical model results 
are from Botvina et al. |12| at a density of 10 13 g/cm 3 . 



nuclei (A) as 300. Average (^4) and (Z) are collected in 
Table |U for these two models. The dispersion in the ratio 
of weak to electromagnetic charge N/Z is important for 
neutrino scattering. This ratio is plotted in Fig. [21 and 
the dispersion is more similar for the two models although 
somewhat larger for the microscopic model. 
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FIG. 2: Ratio of neutron number to charge N/Z versus Z 
for ions at a temperature of T = 1 MeV and proton fraction 
Yp = 0.2. The microscopic Pasta results are from ref. |14| at 
a density of 1.66 x 10 13 g/cm 3 while statistical model results 
are from Botvina et al. |12| at a density of 10 13 g/cm 3 . 



FIG. 3: Distribution of charge Z versus neutron number N 
for ions in a statistical model [T^ at a density of 10 12 g/cm 3 , 
temperature of T = 1 MeV and proton fraction Y p — 0.5. 



B. Static Structure Factor 

In this subsection we present molecular dynamics simu- 
lation results for the static structure factor S(q). We ran 
two simulations with 4000 ions. The first, labeled CHI 
in Table [HI was based on the mixture in Fig. ^ from 
the pasta results of 14]. The second, labeled CH2, was 
for a pure system assuming a single component plasma 
where each ion has neutron number iVj = (N) and charge 
Zi = (Z) equal to the average values of the mixture. Both 
simulations were done at density 1.66 x 10 13 g/cm 3 , pro- 
ton fraction Y p = 0.2, and temperature T = 1 MeV. 
For these conditions the ion density is pi — 7.18 x 10~ D 
fm~ 3 , see Tabled The mixture simulation was started 
from random coordinates and warmed up for a time of 
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FIG. 4: Ratio of neutron number to charge N/Z versus Z 
for ions in a statistical model at a density of 10 12 g/cm 3 , 
temperature of T = 1 MeV and proton fraction Y v — 0.5. 



TABLE II: Simulation runs: the composition is from a mi- 
croscopic model (Pasta) [l4j . a statistical model of Botvina 
(Bot.) 12] or a single (pure) species, the time step is At, while 
Tw is the warm up time and Tm the measurement time, and 
the potential energy per ion is (V) 



Run 


■AT- ion 


Comp. 


At 


p;(10" 5 


7V(10 5 


T M (10 6 


(V) 








fm/c 


fm- 3 ) 


fm/c) 


fm/c) 


(MeV) 


CHI 


4000 


Pasta 


50 


7.18 


5 


1.6 


367.467(1) 


CH2 4000 


Pure 


50 


7.18 


1.6 


1 


368.94(1) 


CH3 


1000 


Bot. 


100 


1.06 


4.1 


8 


162.952(2) 


CH4 


1000 


Pure 


100 


1.06 


2.4 


8 


161.990(3) 


LCI 


1000 


Pasta 


56 


7.18 


2.8 


1.1 


357.65(9) 


LC2 


1000 


Pure 


56 


7.18 


2.8 


1.1 


359.08(8) 


DB1 


40000 


Pasta 


50 


7.18 


1.8 


13 


367.712(1) 


DB2 


10000 Pure 


50 


7.18 


14 


5 


368.151(1) 



500000 fm/c. The system was then evolved for a further 
1.6 x 10 6 fm/c during which S(q) was measured from con- 
figurations written to disk every 500 fm/c, see Table ITT1 
The pure system was also started from random coordi- 
nates, and warmed up for a time of 1.6 x 10 6 fm/c. It was 
then evolved for 1.0 x 10 6 fm/c, with S(q) again measured 
from configurations written out every 500 fm/c. S(q) for 
these runs are compared in Fig. [S] In Fig. E]we compare 
them for small q, as well as with the Debye Huckel results 
ofEq. Eg) . 

The first thing we note from these figures is that S(q) 
for the pure system is systematically smaller than S(q) 
for the mixture. This is because the mixture can have 
additional fluctuations in the weak charge density from 
interchanging isotopes that do not have corresponding 
fluctuations in the charge density, see below. 

Second, Figure El shows good agreement between the 
MD simulation results for S(q) and the Debye Huckel 



FIG. 5: Static structure factor S(q) versus momentum trans- 
fer q for a system at a density of 1.66 x 10 13 g/cm 3 , Y p = 0.2, 
and T — 1 MeV. The solid curve uses a distribution of ions 
from a microscopic model (CHI) while the dashed line as- 
sumes a single component plasma (CH2). Results are from 
molecular dynamics simulations using 4000 ions. Note the log 
scale. 



approximation S® H , Eq. for q between 0.025 and 
0.055 fin" 1 . For q > 0.055 fm" 1 , S° H is smaller than 
S{q) from MD. Therefore, S° H is only valid for small 
q. The lowest q MD result in Fig. corresponds to 
q = q-min = 2n/L where L 3 is the simulation volume. 
There may be a systematic error for this point associated 
with the finite measurement time of 1.6 x 10 6 fm/c for 
the MD calculation of S(q). It may take a long time 
for fluctuations in the weak charge to diffuse across the 
system and this may require a long measurement time to 
obtain an accurate S(q) at low q. We discuss this further 
in subsection IIV Dl 

We also performed MD simulations comparing the sta- 
tistical mixture of Botvina against a pure system. These 
are labeled CH3 and CH4 in Table [HJ and involve just 
1000 ions. Rather using the Botvina distribution shown 
in Fig. ^ however, where p = 10 13 g/cm 3 and Y p = 0.2, 
we used the mixture of Figs. |3| EI where p = 10 12 
g/cm 3 and Y p = 0.5. We kept the temperature at 
T = 1 MeV. The reason is that the denser system may 
be a solid. The ratio of a typical Coulomb to ther- 
mal energy is T = Z 2 a/aT where the ion sphere ra- 
dius is a = (3/47rpi) 1 / 3 . Assuming a pure system with 
Z = (Z) = 75.7 we have T = 335. A pure one component 
plasma is expected to be a solid for T > 180. Note, the 
behavior of a mixture of ions could be somewhat different 
from a pure system. However, given the very large value 
of r the mixture could still solidify. 

The S(q) curves of runs CH3 and CH4 are shown in 
Figs. [7]and|Hl Now there is agreement between MD and 
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FIG. 6: Detail of the static structure factor S(q) versus mo- 
mentum transfer q at low q from Fig. |5] Debye Huckel results, 
Eq. 1161 for a mixture of ions are shown as the solid line while 
the dashed line shows results for a pure system. Note the log 
scale. 



FIG. 7: Static structure factor S(q) versus momentum trans- 
fer q for a system at a density of 10 12 g/cm 3 , Y p = 0.5, and 
T — 1 MeV. The solid curve uses a distribution of ions from a 
statistical model (run CH3) while the dashed line assumes a 
single component plasma (run CH4). Results are from molec- 
ular dynamics simulations using 1000 ions. Note the log scale. 



Debye Huckel results near q « 0.015 fm _1 . For Y p = 0.5 
the spread in N/Z for the ions is smaller than at Y p = 0.2. 
As a result, S(q) for the mixture is enhanced over S(q) for 
the pure system by a smaller amount than at Y p = 0.2. 
These results for S(q) will be used in subsection IIV El 
to calculate neutrino mean free paths. However, first 
we present results for the dynamical response function to 
gain insight into the difference between S(q) for mixtures 
and S(q) for a pure system. 



C. Dynamical Response Function 

In this subsection we calculate the dynamical response 
function S(q, w) to study further the difference between 
pure systems and systems with mixtures of ions. The 
dynamical response function describes the probability for 
a neutrino to transfer momentum q and energy w to the 
system. We have calculated S(q,w) for a microscopic 
nucleon model in ref. 01- The static structure factor is 
the energy integral of S(q,w), 

/>oo 

S(q) = / S(q,w)dw. (25) 
Jo 

We calculate S(q, w) as an integral over the density- 
density correlation function S(q, t) as follows, 

1 [ Tma * 

S(q,w) = - S{q,t)cos{wt)dt, (26) 

Jo 



where S(q, t) is, 

S(q, t) = - 1 / p(q, t + s)*p(q, s)ds. (27) 

+ * ion-*-ave JQ 

Here p(q,t) is p(q), Eq. evaluated with ion coordi- 
nates Ti(t) at time t. We note that S(q) = S(q,t = 0). 
Figure El shows S(q, w) for q = 0.116 frn -1 at a density of 
1.66 x 10 13 g/cm 3 , Y p = 0.2, and T = 1 MeV, assuming 
either a mixture of ions from the microscopic model or a 
single ion species. All simulations show a large peak near 
W = 0.003 frn -1 (0.6 MeV) that corresponds to plasma 
oscillations of the ions 0] . This peak is virtually iden- 
tical for calculations with a mixture of ions or a single 
species. Plasma oscillations may depend on the ratio of 
charge density to average ion mass, since the restoring 
force depends on the charge density while the oscillation 
frequency also depends on one over the ion mass. How- 
ever plasma oscillations do not appear to depend on the 
dispersion in the ion charge or mass distributions. 

In addition, the simulation for a mixture of ions shows 
a large peak at w — 0. This peak is absent in simulations 
with a single ion species. We conclude that this w = 
peak is the primary difference between mixtures and a 
single species. The static structure factor of a mixture 
is larger than S(q) for a pure system by the area under 
this w — peak. Fluctuations in weak charge density, at 
fixed charge density, feel no electrostatic restoring force. 
Therefore these fluctuations will diffuse slowly through- 
out the system and contribute to the response at low 
w. However, such fluctuations in weak charge concentra- 
tion are only possible for mixtures. In a pure system the 
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FIG. 8: Detail of the static structure factor S(q) versus mo- 
mentum transfer q at low q from Fig. Q The circles are 
results for a mixture of ions while the squares are for a pure 
system. Debye Huckel results, Eq. 11611 for a mixture of ions 
are shown as the solid line while the dotted line shows results 
for a pure system. Note the log scale. 



weak charge density must be proportional to the elec- 
tromagnetic charge density and there can be no weak 
charge fluctuations without electrostatic restoring forces. 
Therefore the response of the pure system has no large 
peak at w = 0. 

The peak at w = is also seen in the microscopic calcu- 
lations of Ref . ■ Here molecular dynamics simulations 
were preformed directly in the nucleon coordinates. Nu- 
cleons dynamically formed clusters (nuclei) of different N 
and Z under the influence of nuclear and Coulomb inter- 
actions. Fluctuations in concentration of these clusters 
can then diffuse as in our ion simulations. In addition, 
these microscopic calculations have two effects that are 
not in the ion simulations. First, the nuclei have a variety 
of excited internal modes. Second, nucleons can diffuse 
into and out of clusters leading to nuclear reactions and 
the changing of N and Z. These microscopic degrees of 
freedom may be responsible for significantly broadening 
the plasma oscillation peak. However, since both the ion 
and microscopic simulations have the w — peak, we 
conclude that these microscopic degrees of freedom are 
not necessary for the low energy mode. 

Finally, we have calculated S(q,w) for a pure system 
by approximating the integrals in Eqs. (|26I27|I by a sum 
using a time step At of 250 or 50 fm/c. In both cases the 
ion trajectories have been integrated with a time step of 
50 fm/c. At high w near 0.005 fm the calculation with 
At = 250 fm/c shows a small but nonzero value. This 
appears to be a numerical artifact because the response 
for At = 50 fm/c is much smaller. However, the response 
at w — does not show a peak for either value of At. 
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FIG. 9: Dynamical response function S(q, w) versus energy 
transfer w for q — 0.116 fm -1 at a density of 1.66 x 10 13 
g/cm 3 , Y p = 0.2 and T = 1 MeV, from MD simulations with 
4000 ions. The solid curve assumes a mix of ions from the 
microscopic model fl4|. and has a large peak near w = 0. 
The dotted and dashed curves are for a single component 
pure system using a time step At of either 50 or 250 fm/c 
(see text). 



We conclude that the large w = peak seen in the mix- 
ture simulation, but not seen in either pure simulation, 
is unlikely to be a numerical artifact. 



D. Finite Size Effects 

In this subsection we study the dependence of the MD 
results on the number of ions in the simulation. In figure 
EI]wc compare S(q) results for simulations using 1000, 
4000, and 40000 ions. All simulations are for a mixture 
of ions at a density of 1.66 x 10 13 g/cm 3 , as in Fig. [5] For 
momentum transfers below q = 0.055 fm -1 there is good 
agreement between MD simulation results for 40000 ions 
and the Debye Huckel approximation. However, despite 
measuring for the relatively long time of 1.3 x 10 7 fm/c, 
see run DB1 in Table ITT1 there are still significant statis- 
tical errors at the lowest q values. A MD simulation with 
4000 ions, see run CHI of Table ITT1 agrees with the De- 
bye Huckel S(q) (and with run DB1) for 0.02 < q < 0.055 
fin . However, the smallest q point near 0.017 fm is 
somewhat lower. Finally, a simulation with 1000 ions, 
see run LCI of Table ITT1 underpredicts the Debye Huckel 
S(q) for q < 0.05 for 1 . 

In order to minimize finite size effects, we calculate 
S(q) for q greater than or equal to a minimum value 
q-min — 27r / L with L the box size. This minimum q is 
clearly smaller for simulations with more ions. However, 
the relatively large statistical errors found for S(q) at 
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FIG. 10: Static structure factor S(q) versus momentum trans- 
fer q for conditions as in Fig. |S] assuming a mixture of ions 
from a microscopic model |l4JI . Molecular dynamics results 
for 1000 ions are solid triangles, 4000 ions are open squares, 
and results for 40000 ions are filled circles. The solid curve is 
the Debye Huckel approximation. 



FIG. 11: Static structure factor S(q) versus momentum trans- 
fer q for conditions as in Fig. |S] assuming a single ion species. 
Molecular dynamics results for 1000 ions (triangles) , 4000 ions 
(squares) and 10000 ions (circles) are shown. The solid curve 
is the Debye Huckel approximation. 



E. Neutrino Mean Free Path 



small q may be because it takes a long time for concen- 
tration fluctuations to diffuse across a large box size L. 
Therefore, in order to simulate S(q) accurately at small 
q it appears to be necessary to both use a large number 
of particles and to measure for a long time. Nevertheless, 
we conclude from Fig. 1101 that molecular dynamics sim- 
ulation results for S(q) do converge to the Debye Huckel 
approximation for small q and that this convergence is 
clearly seen in simulations with 4000 or 40000 ions and 
for q < 0.055 £m -1 . 

Figure ITU shows S(q) for a pure system composed of a 
single ion species from MD simulations using 1000, 4000, 
and 10000 ions, see the LC2, CH2, and DB2 runs respec- 
tively in Table [H] Finite size effects appear smaller than 
in Fig. ^3 For a pure system there are no fluctuations 
in concentration. Therefore difficulties with large diffu- 
sion times and statistical errors both appear reduced at 
small q. Again the MD results appear to converge to the 
Debye Huckel approximation, although this convergence 
may only occur at smaller q < 0.03 fm^ 1 compared to 
the q < 0.055 fin" 1 of Fig. UJ 

The agreement between our MD results and the Debye 
Huckel approximation at low q provides a significant test 
of our simulation codes. Furthermore three independent 
MD codes CH (used for runs CHI through CH4 of Table 
llll etc.) , DB, and LC were used. There appears to be 
good agreement between the codes. 



In this subsection we use our S(q) results to calculate 
neutrino scattering mean free paths. Figure lT2*l shows the 
neutrino transport mean free path A (for nucleus elastic 
scattering) versus neutrino energy at a density of 1.66 x 
10 13 g/cm 3 based on simulations CHI, and CH2 of Table 
UTI We see that ion-ion correlations significantly increase 
A compared to that for free ions. 

Next, Fig. 1131 shows the ratio of mean free paths from 
Fig. 1121 for a mixture of ions compared to a single ion 
species. The mean free path for a mixture can be up to 
a factor of two shorter than that for a pure system. 



V. SUMMARY AND CONCLUSIONS 

Neutrinos in core collapse supernovae are likely first 
trapped by large neutrino-nucleus elastic scattering cross 
sections. These cross sections are reduced by ion-ion cor- 
relations. In this paper, we present classical molecular 
dynamics (MD) simulations of ion systems with strong 
Coulomb interactions. We find that neutrino cross sec- 
tions for mixtures of ions, with a dispersion in the ratio of 
neutron to proton number, are systematically larger than 
those for a pure system composed of a single ion species. 
We consider ion compositions from both a microscopic 
dynamical model and a statistical model. 

To investigate the difference between ion mixtures and 
pure systems we calculate the dynamical response func- 
tion S(q, w). This describes the probability for a neutrino 
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FIG. 12: Neutrino transport mean free path A versus neutrino 
energy E v at a density of 1.66 x 10 13 g/cm 3 . The solid curve 
shows results for a mixture of ions while the dashed curve is 
for a pure single species of ion. Finally, the dotted curve is 
for free ions and neglects ion-ion correlations. 



to transfer momentum q and energy w. We find that mix- 
tures have an extra peak in S(q, w) at w = that cor- 
responds to diffusion of composition fluctuations. This 
peak is absent in simulations of a single ion species. 



Our exact MD simulation results for the static struc- 
ture factor S(q) reduce to the Debye Huckel approxima- 
tion in the limit of small momentum transfer q. However, 
this reduction may only happen at very small q. We have 
studied finite size effects by comparing MD simulations 
with 1000 to 40000 ions. Finite size effects appear small 
for 4000 or more ions. However, it may be necessary to 
measure S(q) for long simulation times in order to obtain 
accurate results at small q. 



Neutrino transport mean free paths can be as much as 
a factor of two shorter in mixtures compared to a pure 
system, for low neutrino energies and neutron rich con- 
ditions. The dispersion in composition may be smaller 
in less neutron rich systems. Finally, an important re- 
maining uncertainty is simply the average ion composi- 
tion predicted by different models. The neutrino mean 
free path can be significantly influenced by changes in 
the average nuclear size. 
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FIG. 13: Ratio of neutrino transport mean free paths versus 
neutrino energy E v for a mixture of ions compared to that for 
a pure system. Conditions as in Fig. E2 
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